#!/usr/bin/perl

use strict;
require "CreateCDS.pl";



############################################################################
# Read user inputted restriction sites
#########################################
# Sites must be more than 5 nucleotides long.
# Both fwd and rev versions of restriction sites will be checked against the newly generated sequence.
# If this program fails to generate a sequence with no restriction sites, a warning will be sent to terminal.
# The shorter the site length, the more likely the program will fail (pretty unlikely for sites > 5 nt).
############################################################################

my @rest_sites = getRestrictionSites("sites.fasta");
#print "\nRestriction sites:\n@rest_sites\n";

#########################################################################
# end read restriction sites
#########################################################################








# Riboswitch and Shine-Dalgarno/RBS binding site
##########################################################################

# The effector sequence is a user inputted DNA string that will bind the riboswitch and release the Shine-Dalgarno ribosome binding site from an inhibitory hairpin. The effector sequence must be >= 15 nt. This sequence will be referenced against the user inputted restriction sites to check for conflicts. If a user inputted restriction site is discovered within the user inputted effector sequence, a warning will be printed to terminal.

##########################################################################

my $effector = getSingleFastaSequence("effector.fasta");
$effector = lc $effector;
$effector =~ tr/u/t/;
my $obs = $effector;
$obs =~ tr/acgt/tgca/;
$obs = join("", reverse( ($obs =~ m/./g) ));

# check for restriction sites in effector sequence

my $obs_flag = check_restrict_sites($obs);
if ($obs_flag == 0){
  print "\n *** Warning: Effector sequence contains restriction sites ***\n";
};


# end effector sequence read
#################################
my $iterations_thus_far=0;
my $seed = getSingleFastaSequence("protein.fasta");
my @parameters = getParameters("params.txt");
my %codonFreqs = getCodonFreqs("codonfreq.txt");
open OUTPUT, ">output.fasta" or die "couldn't open output.fasta: $!\n";



START:
# Initialize global variables
my $error=0;
my $obs_error=0;
my $jump = 0;
my $effect_check = 0;
my @nuc_array = qw( a g t c);

#################################

my $obs_len = length $obs; # effector binding site length

=no limit on obs length
if ($obs_len < 15){
  print "\n*** Effector must be 15 or more nucleotides long for riboswitch generation ***\n\n"; 
  die;
};
=cut


# Generate riboswitch
#########################################################
#########################################################
RIBOSWITCH:

my $flag = 0;
my $counter = 0;
my $rbs = 'aggaggt'; # e.coli Shine Dalgarno ribosome binding site - invariable;
my $rbs_len = length $rbs; # length of ribosome binding site
my $rbs_origin = 0; # subtring origin -- changed if restriiction site is found. This origin is important for matching between the RBS complement and the 2nd (upstream of 1st) RBS site



VAR_R:

### Complementary hairpin sites 
##########################################


my $rbs_2 = substr($rbs, $rbs_origin,4);
my $rbs_comp = reverse( substr($rbs, $rbs_origin,4) );
$rbs_comp =~ tr/agtc/tcag/;
my $rbs_comp_length = length $rbs_comp;

my $begin = substr($obs,3+$rbs_origin,2);
my $obs_rev = reverse $effector; # 3'->5' version of effector / complementary effector binding site (obs)
my $post_rbs = substr($obs_rev,(-5-$rbs_origin),5); # complementary effector binding site residue downstream of ribosome binding site
my $post_rbs_comp = substr($obs_rev,(-2),2); # complementary effector binding site residue downstream of ribosome binding site
my $post_len = (length $post_rbs) ; # length of above


my $rbs_obs_len = $post_len + $rbs_comp_length + $rbs_origin  ; # computed spacer on effector binding to skip ribosome binding site (so complementary bases match up even with the rbs space)
my $pre_rbs = reverse( substr( $obs , $rbs_obs_len , 4) ); # effector binding residue upstream of ribosome binding site
$pre_rbs =~ tr/agtc/tcag/; 
#$pre_rbs_comp =~ tr/agtc/tcag/;


## Add spacer after RBS and before ATG 
######################################

my @special_nuc_array = qw( a a g g t t t c);
my @var4 = ();
for( my $i = 0; $i < 7; $i++) {  
  my $rand= join '', map $special_nuc_array[rand @special_nuc_array], 1;
  $var4[$i] = $rand;
};

my $var4 = join($_,@var4); # transform generated DNA sequence to string

if ($var4 =~ m/atg/gi){goto VAR_R};

##################################################
my $riboswitch1 = $begin . " " .$rbs_comp ." " . $obs ." " . $pre_rbs . " " . $rbs_2 ." " . $post_rbs ." " . $rbs ." " . $var4;

my $riboswitch =  $begin . $rbs_comp . $obs . $pre_rbs .  $rbs_2 . $post_rbs . $rbs . $var4;


my $stem_ref = $riboswitch;
$stem_ref =~ s/$obs/"x" x $obs_len /e;
$stem_ref =~ tr/atgc/./;
$flag = check_restrict_sites($riboswitch);
$counter++;

if ($counter > 100){
  $error++;
  $flag = 1;
}
print "\n\n  riboswitch:\n  $riboswitch1\n\n"; #constraint:\n $stem_ref\n\n";
if($flag =~ 0){
  $rbs_origin = int(rand(2));
  goto VAR_R;
}


if ($jump==1){goto CHECKER};
##############################################################################
# end effector / rbs site generation
###########################################################################








#########################################################################
# Promoter Generation
#####################
# Promoter region will be generated with an invariable -35 element and an invariable Pribnow box.
# Variable regions will be generated between segments for proper spacing
##########################################################################

PROMOTER:

$flag = 0;
$counter = 0;
my $obs_counter = 0;
my @var1 = ();
my @var2 = ();
my @var3 = ();
my $promoter = '';
my $var1 = '';
my $var2 = '';
my $var3 = '';
my $box35 = 'ttgaca'; #TTGACA; # consensus -35 promoter element
my $pribnow = 'tataat'; #TATAAT; # consensus -10 pribnow box promoter

while ($flag == 0){


  @nuc_array = qw(a t g c); # array carrying nucleic acids


# variable regions (1, 2, 3) in promoter sequence flanking promoters
###########################################################################

# variable region upstream of -35 element
    for( my $i = 0; $i < 10; $i++) {  
      my $rand= join '', map $nuc_array[rand @nuc_array], 1;
      $var1[$i] = $rand;
    };

# variable region between promoter elements
  for( my $i = 0; $i < 19; $i++) {  
    my $rand= join '', map $nuc_array[rand @nuc_array], 1;
    $var2[$i] = $rand;
  };

# variable region downstream of pribnow box, at start of transcription

  for( my $i = 0; $i < 4; $i++) {  
    my $rand= join '', map $nuc_array[rand @nuc_array], 1;
    $var3[$i] = $rand;
  };

# push variable region arrays into sequence strings
  $var1 = join($_,@var1); 
  $var2 = join($_,@var2); 
  $var3 = join($_,@var3); 

# combine regions into upstream promoter 
  $promoter = $var1 . $box35 . $var2 . $pribnow . $var3;


  $flag = check_restrict_sites($promoter);
  if ($counter > 100){
    $error++;
    last;
  }
  my $effect_check = check_effector_sites(my $terminator_string);
  if ( $effect_check == 0){
    $flag = 0;
    $obs_counter++;
  }

  if ($obs_counter > 100){
    $obs_error++;
    last;
  }  

}

print "\nPromoter:\n $var1 $box35 $var2 $pribnow $var3\n";
if ($jump==1){goto CHECKER};
# end promoter region generation
##########################################################################
##########################################################################





# test protein sequence
# Stephen's protein generator goes here
############################################################################


# variable region downstream of pribnow box and start of transcription
CDS:
$flag = 0;
$counter = 0;
$obs_counter = 0;
my $seq = '';
my $cds;

my $matches_rest;
my $matches_effector;
my $premature_termination;
do {
  $cds = getCodingSequence($seed, \@parameters, \%codonFreqs, \@rest_sites, $effector);
  $matches_rest='false';
  $matches_effector='false';
  $premature_termination='false';

  my @cods = ($cds =~ m/.{3}/g);
  pop @cods;
  for my $cod (@cods) {
    if($cod =~ m/TAA|TAG|TGA/i) {
      # print "premature termination found; regenerating...\n";
      $premature_termination='true';
    }
  }

  for my $rs (@rest_sites) {
    $rs =~ tr/uU/tT/;
    if ($cds =~ m/$rs/i) {
      #print "Sequence matches $rs; regenerating...\n";
      $matches_rest='true';
    }
  }
  if ($cds =~ m/$effector/i) {
    #print "Sequence matches $effector; regenerating...\n";
    $matches_effector='true';
  }
} while (($cds eq 'noseq') or ($matches_rest eq 'true') or ($matches_effector eq 'true') or ($premature_termination eq 'true'));

#print "\nCoding Sequence:\n", $cds, "\n";

my @seq = ();
my $terminator_string = '';

while ($flag == 0){
  for( my $i = 0; $i < 100; $i++) {  
    my $rand= join '', map $nuc_array[rand @nuc_array], 1;
    $seq[$i] = $rand;
  };

  $seq = join($_,@seq); 

  $flag = check_restrict_sites($seq);
  if ($counter > 100){
    $error++;
    last;
  }

  my $effect_check = check_effector_sites($terminator_string);
  if ( $effect_check == 0){
    $flag = 0;
    $obs_counter++;
  }

  if ($obs_counter > 100){
    $obs_error++;
    last;
  }  


}

print "\nCoding Sequence:\n$cds\n";

#end protein coding sequence generation
############################################################################
############################################################################




############################################################################
# Terminator Sequence
# 
# Enters a while statment with a flag set to zero and generates an
# intrinsic terminator sequence. Only exits the while loop when the sequence
# generated does not contain any restriction sites
##
############################################################################
TERMINATOR:
$flag = 0;
$counter = 0;
$obs_counter = 0;
my @var5 = ();
my $var5 = '';

while ($flag == 0){


  for( my $i = 0; $i < 2; $i++) {  
    my $rand= join '', map $nuc_array[rand @nuc_array], 1;
    $var5[$i] = $rand;
  };
  $var5 = join($_,@var5); # transform generated DNA sequence to string


    $terminator_string ='';
  my $stem_seq = generate_stem();
#print "steme_seq = $stem_seq\n";
  my $stem_compliment = $stem_seq;
  $stem_compliment =~ tr/CG/GC/;
  $stem_compliment = reverse($stem_compliment);
#print "compl_seq = $stem_compliment\n";

  my $terminator_loop = generate_loop();
  $terminator_string = $var5 . 'AA' . $stem_seq . $terminator_loop . $stem_compliment . 'TTTTTTTT';
#print "out of loop: $terminator_string";

# ***********    Jacob turned below off after addition of CHECKER section ************************

### Ensures no restriction sites have been generated in the sequnence
#$temp_length = length($seq);
#$end_of_seq = substr($seq, $temp_length - 7, 7);
#$terminator_w_tail_of_seq = $end_of_seq . $terminator_string;
#$flag = check_restrict_sites($terminator_w_tail_of_seq);

# ***************************************************************************************
  $flag = check_restrict_sites($terminator_string);
  $counter++;
  if ($counter > 200){
    $error++;
    last;
  }

  $effect_check = check_effector_sites($terminator_string);
  if ( $effect_check == 0){
    $flag = 0;
    $obs_counter++;
  }

  if ($obs_counter > 100){
    $obs_error++;
    last;
  }  

}
$terminator_string = lc $terminator_string;
 print "\nTerminator Sequence:\n$terminator_string\n";
if ($jump =~1){goto CHECKER};
# end terminator
###########################################################################
###########################################################################


# Subroutines
############################################################################
# sub generate_stem
# 
# Generates a stem of 6 nt long that is entirely C/G content
##
############################################################################
sub generate_stem{

  my $temp_stem = '';
  for (my $i = 0; $i<6; $i++){
    my $random_number =rand();

    if ($random_number > .5){
      $temp_stem = $temp_stem . 'C';
    }
    if ($random_number < .5){
      $temp_stem = $temp_stem . 'G';
    }
  }
  return $temp_stem;
}

############################################################################
# sub generate_loop
# 
# Generates a Hairpin loop of 4 nt that is rich in U/G/C content but without
# forming Watson-Crick-Franklin basepairs
##
############################################################################
sub generate_loop{
  my @temp_loop =();
  my $j = 3;
  for (my $i = 0; $i<2; $i++){
    my $random_number =rand();

    if ($random_number < .25){
      $temp_loop[$i] = 'T';
      $temp_loop[$j] = 'G';
    }
    if ($random_number > .25 && $random_number < .50 ){
      $temp_loop[$i] = 'T';
      $temp_loop[$j] = 'C';
    }
    if ($random_number > .50 && $random_number < .75){
      $temp_loop[$i] = 'C';
      $temp_loop[$j] = 'T';
    }
    if ($random_number > .75){
      $temp_loop[$i] = 'G';
      $temp_loop[$j] = 'T';
    }
    $j--;
  }

  my $terminator_loop = '';
  for (my $i = 0; $i<4; $i++){
    $terminator_loop = $terminator_loop . $temp_loop[$i];
  }

  return $terminator_loop;
}
############################################################################
#end terminator coding sequence generation
############################################################################



############################################################################
# sub check_restrict_sites
# 
# Looks for restriction sites. If found returns a 0 otherwise if not found
# Returns a 1.
##
############################################################################

sub check_restrict_sites {

  my $input_string = $_[0];

#print "in loop '@rest_sites'\n";
#print "\n\ninput: $input_string\n\n";
  foreach my $element(@rest_sites) {
#print "\n\nelement: '$element'\n\n";
    if ($input_string =~ m/$element/gi){
      #print ">>>>>>>>>> a restriction site has been found:$element.\n";
      return 0;
    }
  }
  return 1;
}

#######################################################################
# end restriction check subroutine
#######################################################################


############################################################################
# sub check_effector_sites
# 
# Looks for effector sites outside of OBS riboswitch region. Returns a 0 otherwise if not found
# Returns a 1.
##
############################################################################

sub check_effector_sites {

  my $input_string = $_[0];

  if ($input_string =~ m/$obs/gi){
#print ">>>>>>>>>> a foreign effector site has been found:$element.\n";
    return 0;
  }

  return 1;
}

#######################################################################
# end effector check subroutine
#######################################################################










# Check restriction sites at junctions of generated sequences
##########################################################################
my $check_count = 0;

CHECKER:
$check_count++;

if ($check_count >50){
  $error++;
  goto EXIT_CHECK;
}

$jump = 1;
my $P_R_junction = $promoter . $riboswitch;
my $R_S_junction = $riboswitch . $cds;
my $S_T_junction = $cds . $terminator_string;

my $P_R_flag = check_restrict_sites($P_R_junction);
if ($P_R_flag == 0){
#print"\n Regenerating PROMOTER-- rest site at junction\n*************\n";
  goto PROMOTER;
  $error=0;
};

my $S_T_flag = check_restrict_sites($S_T_junction);
if ($S_T_flag == 0){
#print"\n Regenerating TERMINATOR-- rest site at junction\n*************\n";
  goto TERMINATOR;
};

my $R_S_flag = check_restrict_sites($R_S_junction);
if ($R_S_flag == 0){
#print"\n Regenerating RIBOSWITCH-- rest site at junction\n*************\n";
  goto RIBOSWITCH;
};


EXIT_CHECK:

# end junction checker
############################################################
############################################################




#Output
###############################################################





my $complete = $promoter . $riboswitch . $cds . $terminator_string;
#print" \nComplete Sequence:\n$complete\n\n";

for my $rs (@rest_sites) {
  if ( $complete =~ m/$rs/i ) {
#    print "matches $rs\n";
  }
}

$complete =~ s/(.{1,80})/$1\n/g;
print OUTPUT ">generated sequence ", $iterations_thus_far+1, "\n";
print OUTPUT $complete;

print "generated sequence ", $iterations_thus_far+1, "\n";

$iterations_thus_far++;
while($iterations_thus_far < $parameters[0]){
  goto START;
}



###############################################################




# Error reporting 
#####################################################

if ($error > 0){
  #print "\n*** Warning: Generated sequence contains restriction sites ***\n";
}


if ($obs_error > 0){
  #print "\n*** Warning: Generated sequence contains multiple effector binding sites ***\n";
}


####################################################
#####################################################

exit;
